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We investigate the constraints on the superfluid fraction of an amorphous solid following from 
an upper bound derived by Leggett. In order to accomplish this, we use as input density profiles 
generated for amorphous solids in a variety of different manners including by investigating Gaus- 
sian fluctuations around classical results. These rough estimates suggest that, at least at the level 
of the upper bound, there is not much difference in terms of superfluidity between a glass and a 
crystal characterized by the same Lindemann ratio. Moreover, we perform Path Integral Monte 
Carlo simulations of distinguishable Helium 4 rapidly quenched from the liquid phase to very lower 
temperature, at the density of the freezing transition. We find that the system crystallizes very 
quickly, without any sign of intermediate glassiness. Overall our results suggest that the experimen- 
tal observations of large superfluid fractions in Helium 4 particles after a rapid quench correspond 
to samples evolving far from equilibrium, instead of being in a stable glass phase. Other scenarios 
and comparisons to other results on the super-glass phase are also discussed. 



I. INTRODUCTION 

Recent experiments on solid He 4 by Kim and Chan [1- 
3] raised, among many others, the important question 
of whether disorder can foster the formation of super- 
fluidity in solid samples. Following earlier theoretical 
analyses [4-6], Ritner and Reppy [7, 8] showed that fast 
quenches produce disordered samples with a change in 
the moment of inertia that corresponds to an extremely 
high fraction of superfluid density, on the order of 20%. 
In addition, the role of He 3 impurities [3] suggest that 
disorder must play an important role in the experiments. 
Other studies [9] suggest that the role of disorder is not 
to enhance the superfluid fraction but instead to induce 
non-equilibrium states in the sample that modify the mo- 
ment of inertia as a function of temperature and fre- 
quency. Consequently, in spite of a long series of theoret- 
ical and experimental studies, the relationship between 
disorder and superfluidity in quantum solids is still not 
clear. 

Here we want to focus on one particular proposal that 
was put forward by Boninsegni et al. [6]: the possibility 
of a bulk long-lived metastable glass phase of He 4 . These 
authors performed a Path Integral Monte Carlo (PIMC) 
numerical simulation of Helium 4 at relatively high den- 
sity (p ~ 0.03 A~ 3 ), where the system was very quickly 
quenched from the equilibrium liquid phase at high T to 
a low temperature T = 0.2 K, at which the HCP solid 
phase is stable. They reported the observation of a phase 
which is structurally similar to the liquid, and with a frac- 
tion of superfluid density as high as 60%; this phase was 
observed to last for a large number of Monte Carlo sweeps 
before the system eventually freezes into the equilibrium 
ordered solid. Boninsegni et al. labeled this the "super- 
glass" phase. Actually, the experimental protocols used 
to solidify Helium likely produce very disordered solids, 
possibly glasses. In fact the experiments in [10] showed 
evidences of very slow dynamics, the hallmark of glassy 



behavior. The natural and still open question is why 
freezing in an amorphous density profile should enhance 
superfluidity compared to the crystalline case, which in- 
stead is thought to show zero or very small condensate 
fractions [11, 12]. Superfluidity is related to exchange, 
which is a local process and depends mostly on the local 
neighborhood of a particle. Thus, one might expect, con- 
trary to the findings discussed above, that dense glasses 
should have a fraction of superfluid density comparable 
to the one of crystals at the same particle density. In- 
deed, a theoretical investigation of the superglass phase 
in a simplified (and yet realistic) model of interacting 
bosons found an extremely small condensate fraction in 
the superglass phase [13]. Clearly, the relation between 
disorder and superfluidity deserves further investigation, 
in order to reach a better microscopic understanding of 
superfluidity in amorphous solids and to explain the nu- 
merical and experimental results. 

The main difficulty in the numerical investigation of 
this problem comes from the fact that the glass phase 
(if any) is always expected to be metastable with respect 
to the crystal phase, which is the true equilibrium phase 
of solid Helium. In a classical system, it is reasonably 
straightforward to get properties of a metastable phase 
or a glass, because one can easily simulate the physical 
dynamics of the system by solving Newton's equations 
of motion [14]. In contrast, the real-time dynamics of 
quantum systems is not accessible numerically because 
of the sign problem, and calculating properties involv- 
ing glassy quantum system is problematic. Previous nu- 
merical work of Boninsegni et al. [6] has looked at the 
fraction of superfluid density of a quenched Helium 4 via 
directly calculating it for a system whose PIMC dynam- 
ics slowly equilibrates. More recently, a quantum version 
of the Mode-Coupling Theory of dynamics in glasses has 
been developed and compared with Path Integral Molec- 
ular Dynamics (PIMD) simulations [15], obtaining ac- 
curate informations on the glass transition in quantum 
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hard spheres. However, in this study exchange effects 
were neglected and therefore superfluidity could not be 
investigated. Therefore, for the moment path integral 
simulations are not conclusive. 

Here we approach the problem in a different way. In 
one of the first works on supersolidity, Lcggett showed 
how one can derive an upper bound for the fraction of su- 
pcrfluid density of a generic many-body system in which 
translational invariance is broken, by means of a varia- 
tional computation [16]. The output of Leggett's compu- 
tation is a formula that needs as only input the average 
density profile of the solid. This formula has been ap- 
plied to Helium crystals, and the aim of this work is to 
use it to study the amorphous solid. At present, there is 
not yet any reliable first principle computation or experi- 
mental measurements of the density profile of amorphous 
Helium 4. We endeavor to generate robust estimates of 
it using a number of different techniques, in particular by 
investigating a model of zero-point Gaussian fluctuations 
around classical configurations, and PIMC simulations 
without exchange (which should be closer to the classical 
dynamics). Checking whether these techniques all give 
roughly similar orders for the bound is a way to assess 
the robustness of our result. In the following, we will 
denote the fraction of superfluid density by "superfluid 
fraction" and we always refer to Leggett's upper bound 
to this quantity, unless otherwise specified. 

The rest of this paper is organized as follows. In sec- 
tion II, we discuss how to adapt Leggett's bound to an 
amorphous solid. In section III A, we compute the bound 
for a profile made of Gaussian fluctuations around a clas- 
sical configuration, and compare the results for an amor- 
phous and an ordered solid, while in section III B we dis- 
cuss previous numerical computations [6]. In section IV 
we try to obtain more precise information by comparing 
a classical simulation of a glass-forming system with a 
PIMC numerical simulation of Helium. In section V, we 
show that under some approximations one can obtain a 
formula for the bound that can - at least in principle - 
be computed from neutron or X-ray scattering data. 



II. LEGGETT'S BOUND 

Leggett showed in his pioneering work on supersolidity 
that the wavefunction of the ground state of a system of 
bosonic particles inside a rotating cylindrical container 
can be obtained by finding the ground state for the non- 
rotating system but with new boundary conditions [16]. 
Using cylindrical polar coordinates and assuming that 
the thickness of the cylinder is much smaller than the 
radius R, the new boundary conditions correspond to 
imposing that the wave function gets an extra phase fac- 
tor exp(— 2itimR 2 ijj /K) when the angle Oi of any particle 
i is shifted by 2ir. Here m is the particle mass and co 
the radial velocity. From the oj dependence of the en- 
ergy of the ground state, E m i„(uj), obtained with these 
new boundary conditions one can compute the superfluid 



density p s by: 

Ps .. 1 d 2 E min (uj) 
— = hm — 

p w->o 1 dw l 

where p is the particle density and Iq = NmR 2 the 
classical moment of inertia. From this expression it is 
clear that upper bounds on the superfluid density can be 
obtained by using variational wavefunctions that in the 
lu — > limit tend to the wavefunction for a non-rotating 
container. Lcggett used a variational wavefunction of 
the form *(ri, •• • ,r N ) = *o(n,--- , f N ) exp[i £\ ip(n)], 
where is the ground state wavefunction for the non- 
rotating case and <f> = J2i a sum °f phases satisfying 
the condition ip(9) = <p(6+2ir)-2TrmR 2 oj/h [16, 17]. The 
bound can be improved by including two-body correla- 
tions [18]. Defining 

P (r}= f dn---df N \Mn,--- ,rW)| 2 ]T<S(f-r;), (1) 
J i 

which is the density profile in the ground state, one finds 
that the variational estimation of E min (oj) reads: 

E min (u) = E o + ^j dr[V V (r)} 2 p(f), (2) 

where Eq is the ground state energy in the non-rotating 
case. 

Because of the assumption that the thickness of the 
cylinder is much smaller than the radius, one can simplify 
the problem even further by "unrolling" the annulus and 
consider the system inside a parallelepiped of length L = 
2ttR in the x direction. In this geometry the phase <p has 
to satisfy the boundary condition ip(0 7 y, z) — <p(L, y,z) — 
voL where vq = mRus/h. The minimization of (2) with 
respect to ip leads to the equation for <p(r): 

V-[p(f)V^(r)] = (3) 
and results in an upper bound on the superfluid density: 

Ps = ^iJ v drp(r)\V<p(7)\ 2 . (4) 

Note that if <p Vo {r) is a solution of (3) with bound- 
ary conditions ip(0,y,z) = ip(L,y,z) — vqL, then ip v i — 
(v' /vo)(p Va is a solution with boundary conditions cor- 
responding to v' . Hence, Eq. (4) docs not depend on 
Vo and we can choose Vo = 1 without loss of generality. 
Furthermore, while in the geometry described above the 
wavefunction should satisfy hard wall conditions at the 
boundary of the box in the y and z directions, we will 
simplify the problem by considering periodic boundary 
conditions in the y and z directions [19]. 

In order to find a solution of Eq. (3) satisfying the 
correct boundary condition is useful to rewrite ip as 

f{r) = v Q ■ f + 5ip(r), (5) 
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where 6<p(r) is defined inside the volume V and satisfies 
periodic boundary conditions, and Vq is a unit vector. 
In the original problem vq — x, but since we reformu- 
lated the problem in a periodic cubic box, the direction 
of vq can be varied without affecting the result, in the 
limit V — > oo. Since 5ip(r) is periodic, we can write the 
equations in Fourier space (see Appendix A for details): 

Q ■ V Pq = ^{q-p)Pq-pl5ipp , (6) 

and from the solution for ibtpq one can obtain the Leggett 
bound [17], that reads in Fourier space: 

7 = 1 "^E^' Q)^qP-? ■ (7) 

P P v o 

Given the density profile, the linear equation (6) for iSip^ 
can be solved by truncating the sum over momenta at 
a given cutoff, \q\ < q max , so that the problem reduces 
to solving a finite set of linear equations, which can be 
done by matrix inversion. We accomplish this via a LU 
decomposition [20]. 

An important remark is that the truncation preserves 
the variational nature of the computation. Indeed, it 
can be seen as setting 5ipg — for \q\ > q ma x, which 
amounts to a particular choice of the variational function 
8tp(f) and hence still gives an upper bound on the true 
superfluid fraction. 

Another important remark is that the bound derived 
above applies only, strictly speaking, to the true ground 
state of the system. In the following however, we are in- 
terested in applying it to the glass state, which is at best 
a long-lived metastable state, the crystal being always 
the true ground state. Still, it is clear from the deriva- 
tion that if the life time r of the state is very long, such 
that for any experimentally accessible frequency one has 
lot 1, then the system does not have time to escape 
from the metastable state during the experiment and the 
bound should apply without modification. 

III. SUPERFLUID FRACTION OF 
AMORPHOUS SOLIDS 

A. Hard sphere systems 

In order to understand whether disorder in the density 
profile can lead to an increase of the superfluid density, 
we shall compare the result of the bound for an amor- 
phous glassy profile and the corresponding crystal. The 
only input for our study are the density profiles of the 
amorphous and crystal state. Unfortunately, the former 
is not available for He 4 in realistic conditions. As a con- 
sequence, we decided for a first study to focus on a more 
simple and academic case that can still provide insights 
on the role of disorder. We consider the amorphous and 
crystalline density profiles that one obtains for classical 




FIG. 1: Leggett upper bound for p s /p, for a Gaussian profile 
of width A 1 / 2 around an amorphous jammed configuration 
and in a FCC lattice, as a function of the adimensional pa- 
rameter £ = p 1 / 3 A 1 / 2 (the Lindemann ratio). 



hard spheres. Although this certainly is not a realis- 
tic model of density profiles for He 4 , it allows us to ad- 
dress the role of disorder on p s . Furthermore, a mapping 
from quantum systems at zero temperature and classi- 
cal Brownian systems allows one to find quantum many 
particle models whose ground state wave-function can be 
mapped exactly on (the square of) the probability distri- 
bution of classical hard spheres systems [13]. Thus, the 
results of this section apply directly to those models. 

Classical hard spheres are known to be characterized 
by a high density crystal FCC phase. However, if com- 
pressed fast enough, or due to a small polydispersity, the 
hard spheres freeze in an amorphous glassy state. A typ- 
ical density profile of a very quickly compressed glassy 
state can be obtained by the Lubachevski-Stillinger com- 
pression algorithm [21] (we used the implementation of 
[22]), which is know to be very efficient in producing 
amorphous jammed configurations. The output of the 
algorithm are the positions 1Z = {Ri,--- , Rn} of the 
particles in a random close packed state (at infinite pres- 
sure). The algorithm is deterministic, but different final 
configurations are obtained by starting the compression 
from random initial configurations of points. The com- 
pression runs were performed at very fast rates (we fixed 
the parameter 7 = 0.1, see [22, 23] for details) in order 
to avoid crystallization. 

Furthermore, we will assume that the density profile of 
a typical glassy configuration at finite pressure is the sum 
of Gaussians centered around the amorphous sites, which 
are the output of the previous algorithm. For classical 
systems, this assumption has been tested numerically for 
FCC crystals [24], and has been often used in density 
functional computations of both ordered [25] and amor- 
phous structures [23, 26], giving accurate results. For 
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quantum systems, the Gaussian model has been shown 
to be accurate enough, at least for the purpose of com- 
puting the Leggett's upper bound [27-29]. 

For a given configuration 1Z, the density profile we use 
is defined as 

p(r\K) =Yl^(\r-Ri\) = J ds lA {\f-8\)Y,5{8-R i ) , 

(8) 

where ja(x) = exp(-|x| 2 /(2^4))/(27rA) 3/2 is a normal- 
ized Gaussian of width A, and \f — Ri\ is the distance on 
the periodic box, i.e. it is the distance between r and its 
closest image of Ri. The corresponding Fourier transform 
reads (neglecting terms of order exp(— L 2 /A)): 

p -(7^ e -V/ 2 _L£ e ^. (9) 

i 

In solving Eqs. (6) and (7) we considered amorphous 
configurations of N = 20 and N — 100 particles. All 
the calculations were done with the cut-off set at q m ax = 
20tt/L. We checked that the result does not depend on 
the specific amorphous configuration used by considering 
different amorphous configurations lZ a , a = I,-- - ,Af; 
this is expected since the superfluid density is a macro- 
scopic quantity. The reported results are therefore av- 
eraged over 10 independent configurations. More details 
on the numerics can be found in Appendix A. 

The results are plotted in Figure 1. One can notice 
that, apart from the smallest values of the dimension- 
less parameter, the two curves corresponding to 20 and 
100 particle configurations perfectly agree. The discrep- 
ancy in the region of small £ = p 1 / 3 'A 1 ! 2 is due to the 
approximation brought by the introduction of a cut-off, 
and vanishes in the limit q m ax 3> l/\/A- 

In order to understand to what extent the disorder in- 
fluences the value of the superfluid density, we compare 
the superfluid fraction found in the amorphous system 
to the values obtained through the same calculations in 
the case of a crystal [17, 27-29]. Figure 1 reports the re- 
sults for the average superfluid fraction of the amorphous 
solid just described and those corresponding to the FCC 
lattice (which is the thermodynamically stable one for 
hard spheres) for the Ri, according to the same Gaus- 
sian model (in the latter case our results are consistent 
with previous ones [17, 27-29]). The difference between 
the two is very small, suggesting two conclusions. 

1. Disorder does not influence much the superfluid 
behavior of the system for comparable values of 
p 1 / 3 A 1 / 2 , at least at the level of this variational 
calculation. 

2. The dependence of p s on the density profile is 
mainly through the Lindemann ratio £ = p 1 ^ 3 A 1 ! 2 . 
This conjecture allows us to obtain an estimate of 
the Leggett upper bound for p s in more realistic 
cases as we will do in the next section. 



To conclude this section, we observe that the above re- 
sults allow to obtain a quantitative upper bound for the 
superfluid fraction of a system whose wavefunction is ex- 
actly the Jastrow wavefunction corresponding to classical 
hard spheres. The quantum glassy phase of this system 
has been discussed in [13]. In both the crystal and glassy 
phases, the values of A 1 ! 2 for classical hard spheres do 
not exceed 0.1 (in units of the sphere diameter) [23-25], 
and the same is true for I, since the density is very close 
to 1 (in the same units) in both solid phases. Using the 
results of Fig. 1, we obtain an upper bound p s /p < 0.1%, 
which is consistent with the extremely small values of the 
condensate fraction found in [13]. 

B. Superfluid fraction of amorphous solid Helium 4 

In this section, we attempt an application of our results 
to the more interesting case of disordered solid He 4 , based 
on the observation above, that an estimate of the Linde- 
mann ratio £ = p x l 3 A x l 2 , together with the results of 
Fig. 1, should provide a reasonable estimate of Leggett's 
bound. 

At the end of Ref.[29] it is stated that, by fitting 
the Path Integral Monte Carlo density profile of HCP 
solid He 4 , one obtains a value \f~A = 0.1274 d at p = 
0.0353 A" 3 and VA = 0.1486 d at p = 0.029 A~ 3 . 
Here d is the nearest-neighbor distance for the HCP 
lattice. The number density of the HCP lattice satis- 
fies the relation pd 3 = y/2, hence d = 2 1 / 6 / / o 1 / 3 and 
£ = \J~Ap 1 l 3 = 2 x ^\J~Ajd. In the same reference it is also 
stated that the upper bound computed by using the fitted 
Gaussian density profile coincides with the one obtained 
by using the true PIMC density profile, and corresponds 
respectively to p s /p = 0.06 and 0.22. These values are 
reported in table I. 

We now make the following assumptions: 

1. At least for the purpose of computing Leggett's up- 
per bound, the true density profile can be fitted to 
a Gaussian profile. This is true for the crystal [29] 
and we assume that it remains true for an amor- 
phous solid. 

2. The parameter £ for the amorphous solid is smaller 
than that of the crystal at the same density. This 
can be understood by observing that crystalline 
configurations are better packed than amorphous 
configurations, therefore leaving more room ("free 
volume") for fluctuations. It is true for Jastrow 
wavefunctions [13] (i.e. classical system) and we 
do not find any reason why quantum fluctuations 
should dramatically affect this property. 

Based on these assumptions, the true Leggett's bound 
for the amorphous system should be smaller than the 
same bound for the crystal at the same density. This can 
be estimated using the values of £ reported in [29] and 
reading the corresponding superfluid fraction from Fig. 1 
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P (A- 3 ) 


HCP, Leggett's bound (Ref.[29]) 

t Ps/P 


Glass, Leggett's bound (this work) 

Ps/P 


Glass, QMC (Ref. [6]) 

Ps/P 


0.029 
0.0353 


0.167 0.22 
0.143 0.06 


0.282 
0.127 


0.6 
0.07 



TABLE I: Leggett's bound for He 4 in the HCP crystal state [29] and glassy state. Quantum Monte Carlo results for the glass 
are also reported [6]. 



or using the results obtained in [29] for the HCP crystal. 
These values are reported in table I and are similar. 

We compare the upper bound obtained in this way with 
the values of p s obtained numerically by Boninsegni et 
al. via PIMC [6]. Interestingly, we find that the bound is 
very close to the PIMC numerical result, and in particu- 
lar at the smallest density the bound is violated by the 
PIMC result. This can be due either to the very rough 
approximations involved in our computation, or to the 
fact that the glass is not a really long-lived metastable 
state at this very low density. The latter possibility, i.e. 
that the system is rapidly evolving out of equilibrium, 
would invalidate the derivation of Leggett's bound but 
it would also raise problematic questions regarding the 
measurement of p s using the Ceperley formula, which is 
strictly valid if thermodynamic equilibrium is achieved 
and in the limit of small frequency. 



IV. DOES A STABLE GLASS STATE EXISTS 
FOR HELIUM 4? 

In order to study the stability of the glass phase in He- 
lium 4, we performed Path Integral Monte Carlo simula- 
tions, that we discuss in this section. Before discussing 
the more complex quantum simulation, we present some 
classical simulations in order to deal with a well con- 
trolled situation, where the presence of a glass transition 
has been firmly established. 



1. What should we expect from a glass-forming system? A 
classical simulation 

We performed standard Molecular Dynamics (MD) 
simulations of the Kob- Andersen binary mixture [14], 
which is known to be a good glass former and does not 
show any sign of crystallization even after very long MD 
runs at low temperature. The latter is a mixture of two 
types of particles (A and B), interacting through different 
Lennard- Jones potentials, with the parameters specified 
in [14] . In the rest of this section we use reduced Lennard- 
Jones units, namely we use uaa and eaa as units of 
length and energy, and m as unit of mass. Consequently, 
yJma' AA /eAA is the unit of time (the latter convention 
is slightly different from the one of [14]). Note that to 
compare with Helium one should keep in mind that for 
that system a <~ 2.56 A and e ~ 10.2 K. 



We quenched a dense (p — 1.2) system of N = 216 
particles from very high temperature (T = 2) to very 
low temperature (T = 0.05) deep in the glass phase (the 
glass transition temperature being around T = 0.435 
at this density [14]). We run the simulation for a to- 
tal time t = 15000 and we printed configurations every 
At = 5 which is of the order of the decorrelation time in 
the glass (estimated from the decay of the self scattering 
functions). From each configuration we deduced 

^) = ^E ei ^' (t) ' ( 10 ) 

3 

where fj(t) is the position of particle j at time t, and the 
corresponding instantaneous value of the static structure 
factor Srft) = V\p^t)\ 2 /p. 

In Fig. 2 we plotted p^(t) and the structure factor S$(t) 
as a function of MD time after the quench. The vectors 
q = 2ir/L(n x , n y , n z ) and the corresponding integers are 
given in the caption. We see that after a short tran- 
sient, the density profiles fluctuate around a non-zero 
value which is quite stable, except for some rare "crack" 
events where the density changes abruptly. These are 
probably due to groups of particles that switch back and 
forth between two different locally stable configurations. 
This system is indeed extremely dense and at very low 
T, therefore its dynamics is basically that of harmonic 
vibrations around local minima of the potential (except 
for the rare cracks). The largest instantaneous value of 
S$(t) corresponds to the (2, 1, —6) curve in Fig. 2 for all 
t > 1000; therefore, all values are smaller than 20 at all 
times, showing that there are no Bragg peaks. This is 
what we expect to see in a glass. In this case, we can 
easily deduce the average values of p^ for a given glassy 
configurations by taking the average of p${t) over a time 
interval where there are no crack events. From these, 
we could compute the Leggett bound as previously dis- 
cussed. 



2. Absence of a stable glass phase from a Path Integral 
Monte Carlo simulation 

Motivated by results of [6] we tried to compute the 
supcrfluid fraction based directly on Path Integral Monte 
Carlo data. Unfortunately, PIMC does not give access to 
the real time dynamics of the system, but following [6] 
we studied the Monte Carlo dynamics, in the hope that 
this is a reasonable proxy for the real time dynamics. 



6 




5000 t 10000 



15000 




15000 




40 



30 



S(q) 

20 
10 



*M*i 



ijiiihiift 



FIG. 2: Evolution of the density profile after a quench from 
high to low temperature for a classical glass forming sys- 
tem, using Molecular Dynamics. (Top) Instantaneous value 
of Sq(t) for three representative values of q (the corresponding 
(n x ,n y ,n z ) are indicated in the caption). (Middle) Instanta- 
neous values of pq(t) for a representative value of q. (Bottom) 
The time average of Sq(t) over the whole simulation, as a func- 
tion of q (in reduced LJ units). Scatter points are values for 
a given q, the full black line is the angular average over all 
vectors with the same modulus. 



FIG. 3: Evolution of the density profile after a quench from 
high to low temperature for a quantum Helium 4 system, 
using Path Integral Monte Carlo. Time here represents the 
number of Monte Carlo sweeps. The panels are the same as 
in Fig. 2, except that the average of S$(t) in the lower panel 
has been taken for t > 75000, and the angular average is not 
reported because of the strong anisotropy of the result. All 
quantities are plotted using A as units of length. 
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The representation of quantum systems in PIMC in- 
volves certain important extensions beyond the classical 
representation of point particles. To begin with, particles 
are represented by paths (or polymers) in space. These 
paths manifest the zero point motion inherent in the 
quantum mechanical system. For distinguishable parti- 
cles, this is the only difference. For particles with statis- 
tics (bosons), these paths then can permute onto each 
other forming larger paths or cycles. 

We initially focus on studying a quenched quantum 
system of Helium particles but require that they act like 
distinguishable particles. There are a number of poten- 
tial advantages of this approach. To begin with, one 
may hope that distinguishable particles are more likely 
to retain the relationship between real dynamics and 
the Monte Carlo dynamics. Secondly, the simulation of 
distinguishable particles is faster and more easily paral- 
lelized over many processors allowing for longer simula- 
tions. 

We used the Aziz potential as a model for Helium [30] , 
and in this section we always use Angstroms as units of 
length and Kelvins as units of temperature. The pair 
product action is used as the approximation for the high 
temperature density matrix and an imaginary time step 
of St = 0.025 K is used. We equilibrated a system of 
N = 216 particles in the liquid phase at a density of 
0.029 A -1 and a temperature of T = 2 K. The system 
is then instantaneously quenched to T = 0.166 K. This 
is accomplished by taking a snapshot of the paths from 
T = 2 K and then, for each time slice of the old path, 
placing 12 time slices for the new lower temperature path; 
this is similar to what was done by Boninsegni et al. [6]. 
We then run the PIMC from this quenched configura- 
tion. These paths are obviously highly artificial because 
the distances between many adjacent time slices are zero. 
Over a very short period at the beginning of the quenched 
run, though, this artificial aspect of the path quickly re- 
laxes leaving the paths in a configuration that mirrors 
the higher temperature formation. 

In the following we refer to t as the PIMC "time" (num- 
ber of PIMC sweeps 1 ), while r is the imaginary time. At 
each "time" t, the PIMC code returns a configuration 
rj(t), the latter being the imaginary time trajectory of 
particle j as function of the imaginary time r. We can 
define the instantaneous density as 

^) = ^Ef dreifrT<<) ' ( n ) 

and the instantaneous structure factor 

S &) = W^f dreW®-™ . (12) 



1 We define a sweep as attempting a displace move on (an ex- 
pected) 10% of the particles and attempting bisection moves on 
(not necessarily unique) 0.47V/(T 8t) time slices. 



Note that in the quantum case, at variance with the clas- 
sical case, these two quantities arc not directly related. 
At each PIMC sweep we recorded the values of the above 
quantities, which we then averaged over 50 PIMC sweeps 
in order to eliminate part of the fluctuations. 

The results for a representative run of the above proce- 
dure are reported in Fig. 3. Unfortunately, the dynamics 
of this system looks quite different from the formation of 
a glass from a quenched liquid. First of all, the structure 
factor becomes quite large for some values of q, therefore 
suggesting the presence of large crystallites in the sam- 
ple. Indeed, the largest value of the structure factor cor- 
responds to the (5, 0, 4) curve in Fig. 3 at large times and 
to the (5, 4, 2) curve in Fig. 3 at short times. We see that 
while at short times the values of S$(t) are smaller than 
10, at larger times they grow up to 50, which clearly indi- 
cates the presence of large crystallites in the sample (note 
in addition that these values have been averaged over 50 
PIMC sweeps and also over imaginary time). Moreover, 
the p$(t) (reported for a representative value of q in the 
middle panel of Fig. 3) are not fluctuating around some 
stable value; they display a sluggish evolution that does 
not allow us to identify a region of times where the system 
is close to some metastable density profile that does not 
evolve in time. What we can learn from this is that the 
quenching from a (exchange-free) liquid to a (exchange- 
free) low temperature liquid froze to a (possibly very 
broken) crystal relatively quickly without showing any 
intermediate signs of glassiness. Note however that this 
behavior was not observed in all runs: some runs did not 
display signs of crystallization for times up to ~ 200000 
PIMC sweeps. Still the dynamics was sluggish enough 
to prevent the identification of a stable glass phase. We 
also tried turning off some moves (the displace moves) in 
order to slow down the relaxation to the crystal, but the 
system still seems to freeze just as quickly. 

In conclusions, we were not able to find a long-lived 
metastable glassy state in our quantum simulations. This 
is probably due to the fact that monodisperse systems al- 
ways crystallize quite fast. This is well known in the clas- 
sical case and seems to also hold true when quantum zero 
point motion is introduced (at least in this specific ex- 
ample) . This leaves the discrepancy between our findings 
and those of [6] to be explained. One possibility is that 
exchange, that we neglected, may be critically important 
for exhibiting the glassy behavior of Helium 4: it could 
be that the path integral at the low temperatures we are 
focusing on is dominated by exchange paths, whereas the 
paths that make the glass unstable are mainly without 
exchange; indeed we find them with our PIMC. In this 
case, the instability of the glass would be a much rarer 
process once one takes into account exchange paths. In 
particular, since crystals have a very low or zero super- 
fluid fraction, we know that their corresponding path in- 
tegral is dominated by paths without exchange. In conse- 
quence, eliminating the exchange could also make crystal 
nucleation easier since it makes it a less rare process. 

An additional possibility is that the glassy behavior is 
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sensitive to the specific details of the simulation (type of 
Monte Carlo moves, length of the paths, etc.). We leave a 
more detailed investigation of this point for future study. 



V. TOWARDS A METHOD FOR 
EXPERIMENTALLY ASSESSING THE LEGGETT 
BOUND 

As we discussed previously, the problem in applying 
our analysis to realistic system is that the amorphous 
density profile of He 4 cannot be easily measured experi- 
mentally. Below, we endeavor to connect the bound on 
p s to the so-called non-ergodic factor g q , which in princi- 
ple could be measured in experiments, e.g. by neutrons 
or X-ray scattering. It is defined as 



N 9q ~ N 



1 \ ~* a a 



(13) 



where the overbar denotes the statistical average over 
the amorphous states sampled statistically by the sys- 
tem. These are indexed by a = 1, • ■ ■ , Af, and under the 
Gaussian approximation each profile p'i is obtained from 
Eq. (9) by plugging the reference positions correspond- 
ing to each different amorphous configuration lZ a . The 
statistical average is performed with the weights a that 
correspond to the frequency with which they appear in 
an experiment, or equivalently their Boltzmann weight. 

First, let us focus on p 8 , which is the average of the 
superfluid density corresponding to each amorphous 
state. Since the superfluid density is a macroscopic quan- 
tity we expect (and we have checked numerically, see Ap- 
pendix A) a self-averaging behavior, i.e. the fluctuations 
of p™ are negligible. However, as usual for disordered 
systems, the computations are easier for ~p s . Multiplying 
Eq. (6) by p°_ - and averaging over a we obtain 



p^Q 



(14) 



where we define, for p,q^0 (that are the only cases 
involved in the equation above) 



Ftf, V) = Tf^2 P^-p iSl PpP-q = Pq-pi^pP-q ■ (1 



5) 



Clearly iipg- is strongly correlated to p^, being the solution 
of (6). In order to simplify the problem we assume that 
these variable are Gaussian distributed. Using Wick's 
theorem, one has 



F{q,T>) = Pq-p iSiPp P-q+ Pq-p iSlffP-^ 



(16) 



+ Pq-piSfp P-q + Pq-pP-q iSfp 

Note that, due to translation invariance of the averages 



over a, one has pq* = p6~$ and Pq-p-p = jj9q5q,p- Hence, 




0.2 1/3 . 1/2 0.3 
P A 



FIG. 4: Result for p s /p as a function of £ = p 1/3 A 1/2 , where 
Ri are the center of the spheres in an amorphous jammed 
configuration of iV spheres with periodic boundary conditions. 
We report the exact computation according to Eq. (6) and the 
approximate result Eq. (20). 



for p,q^0, we get 



F(q,p) = p$-p ibippp-q = pSftfiSipfP-f = pS^pF(q) . 

(17) 

Substituting the last expression in (14), we obtain 



F(q) 



p(q- v )g q 



Nq 2 



(18) 



Averaging (7) over a, we get 



'''' ' VF(q) =1-^2^ „,2„ 2 9q ■ 



q=£0 



v 2 q 2 



(19) 

In the thermodynamic limit, the sum can be replaced by 
an integral, and performing the angular integration we 
obtain: 



P 



dq q 2 
W?P 



9q 



(20) 



The same result can be obtained by means of a large A 
expansion of the system of equations, which however is 
poorly convergent and cannot be used in a systematic 
way, see Appendix B. 

As before we need to introduce a cut-off in the sum on 
q in (9) and calculate numerically the non-ergodic factor 
g q by averaging the density over the same configurations 
lZ a considered above. We set the cutoff according to the 
spherical constraint \q] < q ma x ■ We increased q ma x until 
Imax = 207r/L, when the convergence in g q was reached. 
For the purpose of computing the non-ergodic factor and 
then the approximate bound, as given in Eq. (19), we 
averaged over 100 different configurations. In this case, 
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in fact, one does not face the computational problem of 
inverting the linear system (6) and thus a larger statistics 
can easily be taken. The results of the computations are 
shown in Figure 4. We plotted the superfluid fraction 
obtained through the exact procedure (7) and the ap- 
proximated one (20), both for the configurations with 20 
and 100 particles. The agreement between the approx- 
imated curve and the exact one is good for large value 
of t while they start to differ when the localization pa- 
rameter decreases, for values of the bound around 0.7. 
Unfortunately for the interesting values of i the approxi- 
mated calculation gives wrong results. However, we find 
it useful, since it allows to estimate the typical scale of 
I at which the bound starts decreasing fast from 1 to 
and we hope that it will be possible to improve it in the 
future, in order to be able to apply it to realistic cases. 



VI. CONCLUSIONS 

The aim of this paper was to study Leggett's upper 
bound for amorphous quantum solids. We showed that 
for quantum systems described by a hard sphere Jastrow 
wavefunction, the superfluid fraction must be smaller 
that 0.1%, which is consistent with a previous investi- 
gation that found extremely small condensate fractions 
for this system [13]. Moreover, the hard sphere result 
suggests that crystal and glass phase characterized by 
the same Lindemann ratio should have similar Leggett's 
upper bounds for the superfluid fraction. 

On this basis, we attempted to apply our results to 
glassy He 4 [6]. We found that the upper bound for p s is 
in general very close to the numerical results of Rcf. [6] , 
and at density p = 0.029 A -3 it is below. One possi- 
ble origin of this discrepancy could be that at such low 
density the life time of the metastable glassy state is 
too short, and the system is intrinsically out of equilib- 
rium; in that situation Leggett's bound is inapplicable, 
since it assumes that the reference wave-function corre- 
sponds to a truly metastable state. Indeed we generically 
found from Path Integral Monte Carlo calculations that 
(at least if exchange is neglected) the system crystallizes 
very fast after the quench, which is consistent with a very 
short lifetime of the metastable glass. 

Overall, our findings suggest two possible scenarios 
(not necessarly antithetic). (1) An amorphous stable 
glass has a superfluid fraction, not only a Leggett's upper 
bound, very similar to a defect-free crystal with the same 
Lindemann ratio. Since we know from experiments and 
simulations that this superfluid fraction is very small, 
or possibly zero, we are bound to conclude that the 
glassy supersolid phase found in experiments do not cor- 
respond to a truly stable glass: the system is instead 
rapidly evolving out of equilibrium and, somehow, this 
enhances superfluidity. (2) Exchange promotes glassiness 
and whereas a stable glass phase cannot exist, because it 
has a very short life-time, a superglass can. This could 
be partially tested by comparing the stability of the glass 



phase in imaginary time simulations with and without 
exchange. 

It is worth to note that we neglected the role of small 
concentration of He 3 impurities (of the order of few ppm) 
that has received a lot of attention in experiments [3]. 
The reason is that we focused on a bulk glass phase of 
He 4 , whose density profile should be largely independent 
of such a small concentration of He 3 impurities. It could 
be, however, that He 3 impurities affect the dynamical 
stability of the glass phase. Based on the experience 
on classical systems, it is likely that in presence of a 
large concentration of impurities crystallization will be 
avoided [14] and a long-lived quantum glass phase [15, 31] 
will be stable. In this case, it should be very easy to mea- 
sure the density profile and compute the Leggett bound 
using the procedure detailed above. However, it has been 
estimated that a concentration of at least 0.1% of im- 
purities is needed to stabilize the glass [32]. Therefore, 
the typical concentration of He 3 (<~ ppm) should not be 
enough to produce a sensible effect, unless some unex- 
pected phenomenon related to the quantum mechanical 
nature of the systems (e.g. exchange, as already dis- 
cussed) becomes relevant. 
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Appendix A: Details on the numerical procedure 

We define the Fourier transforms in the cubic box of 
side L and volume V — L 3 as follows: 

P?=XrJ drp{f)e^ , p{r) = £ p q ~e'^ , (Al) 

where q = ^{n x , n y ,n z ), and each of the integers rij G Z, 
and similarly 

5<p 9 = i / drS^(r)e^ F . (A2) 
v Jv 

Note that 5<pq is an irrelevant constant phase in the vari- 
ational wavefunction so we set it to zero. Finally, 

v,= \ V \ *=?' ff (A3) 
which leads immediately to Eq. (6). 
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We performed the calculations for different values of 
the Lindemann parameter £ = p 1 / 3 ^! 1 / 2 , increasing the 
number of vectors q according to the spherical constraint 
\q\ < qmax, until a reasonable convergence in the value 
of the bound (7) was achieved, at least for large val- 
ues of A. From Eq. (9) one sees that for large \q\ the 
corresponding component pq is suppressed through the 

factor e~ Aq I" 1 . Thus, one needs to truncate the sum 
over q at q m ax ~ 1/ VA, as higher terms will not con- 
tribute. Unfortunately, for small A, this cut-off is too 
heavy in terms of computational time and we should use 
a lower one. Still, considering small configurations and 
sufficiently large values of A, which nevertheless span the 
physical region of interest, we could reach a good conver- 
gence or keep the error under control. Note additionally 
that by increasing the number of vectors q in (9), the 
value found for the superfluid fraction monotonically de- 
creases, as expected because of the variational property 
already discussed. This permits to preserve the nature of 
upper bound for Eq. (7), despite the cut-off approxima- 
tion. Overall, we found that the better compromise was 
to set q max = 20-7T/L. 

In order to check the independence of the bound on 
the flow direction, we also compared the results obtained 
with the velocity vo along the (1, 0, 0) direction to those 
along (1,1,1) and we observed a negligible difference 
which is expected to vanish in the thermodynamic limit, 
because amorphous solids are statistically homogeneous 
on large scales. 

We have also checked that the bound for the superfluid 
density almost does not fluctuate by considering different 
amorphous configurations lZ a , a = 1, • • • ,J\f, as it is ex- 
pected since the superfluid density is a macroscopic quan- 
tity. We computed the corresponding superfluid fraction 
p° and the average p s = ^2 a pf/ftf for 10 different config- 
urations. The variance of p s is very small. In this paper 
we presented results averaged over 10 realizations of TZ a , 
a larger statistics do not lead to appreciable differences. 

Finally, as a check of our codes, we repeated all the 
calculations on configurations of 20 particles occupying 
uncorrelated uniformly random positions in the box, i.e. 
where Ri are uniform and independent random vari- 
ables in [0,-L] 3 . In this case it is easy to show that 
g q = cxp(— Aq 2 ). Hence Eq. (20) becomes 

— 2 r°° ^ —Aq 2 1 

7 = 1 ~ 3(2tt)V J qq 6 =1 ~ 24tt 3 / 2 ^ 3 / 2 ' 

(A4) 



In this case the values of the bound were more sensitive 
to the particular realization, so we took averages over 30 
configurations. For every value of the localization pa- 
rameter, the superfluid fractions that we found were on 
average smaller, as reported in Figure 5. 
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FIG. 5: Result for p s /p as a function of the localization pa- 
rameter p 1/la A 1 / 2 , where Ri are N random points in [0, L] 3 
with periodic boundary conditions. 



Appendix B: Large A expansion 

For large A, we expect that the density becomes uni- 
form. Hence, p$ —> p, and pq- — > for q ^ 0. We can use 
this to expand iSipq* systematically in powers of p^. We 
rewrite Eq. (6) as 

q-voPq = q 2 piSLp^+ ^ (q • p)p$-pi8cpp . (Bl) 

We write Stpq* — Stp^ + Sip 1 "? + • • • where the different 
terms are of order (pq) k . At first order 

IP 



at second order 



™V\ ~~— n X^(l' P)PQ-P lSi Pp = ~ 2^1 r,2„2„2 Pf-pPp > ( B3 ) 

at third order 



p^0,q p¥=o,q 



■ (3) 1 /-. • (2) \ - \ - (q ■ P)(p ■ p'W ■ v ) 

% n = —jf~p 2^ ^-p)pt-^y = 1^ — g 2p2 P /2 p3 — pq-pPp-p-pp> (B4) 
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from which we can guess the order k: 

■ _i\fc-i (q-pi)(pi • P2) • • • {pk-i -yp) ,-o-s 

ll fq —\ L > 2-^1 ^2 12 Zk Pq-PlPpi~P2 ' " Ppk-2-Pk-lPpk-l K^^) 

_ , s _ _ , s _ _ _ IPl'-'Pk-lP 

and so on. Plugging this in Eq. (7) we get 

p s _ 1 \ - (vq ■ q) 2 \ " (^0 -p)(p- ^0) 

\ - \ - \ - (go ■ q)(q-p)(p-p')(f ■ go) 

q 2 p 2 p l2 v 2 p 4 Pq-pPp-p'Pp'P-q^ ■ 

While this expansion seems a simple strategy of solution of Eq. (6), it is very poorly convergent and in practice it is 
not very helpful. ^ 
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